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00 ' Abstract. Nacre, commonly known as mother-of-pearl, is a remarkable biomineral 

C^l ■ that in red abalone consists of layers of 400-nm thick aragonite crystalline tablets 

confined by organic matrix sheets, with the (001) crystal axes of the aragonite tablets 
oriented to within ±12° from the normal to the layer planes. Recent experiments 
demonstrate that this orientational order develops over a distance of tens of layers 
from the prismatic boundary at which nacre formation begins. 

Our previous simulations of a model in which the order develops because of 
differential tablet growth rates (oriented tablets growing faster than misoriented 
ones) yield patterns of tablets that agree qualitatively and quantitatively with the 
experimental measurements. This paper presents an analytical treatment of this model, 
O \ focusing on how the dynamical development and eventual degree of order depend on 

model parameters. Dynamical equations for the probability distributions governing 
tablet orientations are introduced whose form can be determined from symmetry 
considerations and for which substantial analytic progress can be made. Numerical 
simulations are performed to relate the parameters used in the analytic theory to 
OO [ those in the microscopic growth model. The analytic theory demonstrates that the 

dynamical mechanism is able to achieve a much higher degree of order than naive 
estimates would indicate. 
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1. Introduction 

Nacre, or mother-of-pearl, is a biomineral that attracts the attention of materials 
scientists, biologists, and mineralogists as well as physicists because of its remarkable 
mechanical properties and its incompletely elucidated formation mechanisms [1-5]. 
Aragonite, a hard but brittle orthorhombic CaCO^ polymorph, accounts for 95% of 
nacre's mass, yet nacre is 3000 times tougher than aragonite [6]. No synthetic composites 
outperform their components by such large factors. It is therefore of great interest to 
understand the mechanisms governing nacre formation. 

Nacre is a layered composite in which organic matrix (OM) sheets alternate with 
aragonite layers, each of which consists of tablets of irregular polygonal shape that 
completely fill the space between preformed OM layers. In red abalone, the OM layers 
are 30-nm thick, and the aragonite tablets are of thickness 400-500nm and width 5-6 
microns [7,8]. The aragonite tablets are crystalline and oriented with their (001) crystal 
axes within ±12° from the normal to the layer plane [8,9] . TEM and AFM measurements 
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have shown that the OM has pores through which aragonite can grow [10, 11], and 
x-ray and x-ray photoelectron emission spectromicroscopy (X-PEEM) measurements 
demonstrate that nacre has stacks of co-oriented tablets [7,8,12], consistent with the 
hypothesis that aragonite crystals grow through pores in the OM sheets. 

Growth through pores explains how crystal orientation of the aragonite tablets 
is maintained between layers in the material, but it does not necessarily explain the 
physical mechanism giving rise to the orientational alignment in the first place. It is 
common belief in the biomineralization field that the alignment of the aragonite crystal c- 
axes is due to microscopic chemical templation by the OM [1,13], with organic molecules 
providing surfaces that promote aragonite nucleation with preferred orientations. 
However, it is not obvious why such a mechanism would lead to a high degree of 
orientational order of the c-axes but not of the a or b axes of the aragonite tablets. 
Moreover, Ref. [8] reports x-ray photoelectron emission spectromicrosopy (X-PEEM) 
and microbeam X-ray diffraction that probe the degree of orientational alignment of the 
aragonite c-axes [7,14] and demonstrate that the orientational alignment of the aragonite 
tablets increases systematically over a length scale of tens of microns starting from the 
prismatic boundary at which nacre growth originates. As discussed in Ref. [8], this 
observation suggests that the aragonite crystal orientation is the result of a dynamical 
process. 

In Refs. [8, 15] a dynamical model is proposed for the development of orientational 
order of the aragonite tablets in nacre in which the presence of co-oriented stacks of 
tablets plays an essential role, and the ordering arises because oriented tablets grow 
faster than misoriented ones. Refs. [8, 15] present numerical simulations of the model 
that successfully reproduce several different aspects of the pattern of tablet orientations 
in red abalone nacre measured using X-PEEM. This paper presents a more detailed 
investigation of the model and its behavior, including closed-form analytic results 
that apply in the limit that the probability of nucleating a tablet with orientation 
different from the one immediately below is small, which, based on comparison to the 
X-PEEM results, is the physically relevant regime. We find that the dynamical ordering 
mechanism can lead to a remarkably high degree of orientation of the tablet c-axes, 
because the width of the distribution of tablet orientations decreases very strongly as 
the fraction of tablets that nucleate with the "wrong" orientation decreases. 

The analytic characterization of the model uses methods similar to those used to 
study mutation-selection models in population biology [16-18]. However, the analytic 
formulation involves some parameters whose relationship to those of the growth model 
is not determined in Ref. [8] . These parameters are examined here, and it is shown that 
this relationship is not trivial. The relationship between the parameter sets depends on 
the spatial arrangement of the nucleation sites, and some aspects are quite insensitive to 
changes in the values of the microscopic parameters. Numerical simulations and mean- 
field arguments are used to relate the parameters in the analytic probabilistic model to 
the parameters in the original growth model for the case when the nucleation sites are 
chosen randomly with uniform probability on each layer. 
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The paper is organized as follows. Section [2] presents the model, while section [3] 
presents the analysis of the probabilistic model that enables one to understand 
qualitatively some features of the behavior, using methods that have been developed 
to study models relevant to population biology [16-22]. The fixed point behavior is 
discussed in subsection 13.11 while the dynamical evolution is discussed in subsection 13.21 
Section H] discusses the relationships between the parameters of the growth model to 
those of the model for the probability distribution for tablet orientations. Sec. [5] is a 
discussion, and the conclusions are presented in Sec. El Appendix A presents additional 
details of the arguments justifying the functional forms of the equations used in the 
main text. 

2. The model 

Figure 1 illustrates the basic mechanism of nacre growth [1,10,23]. First, organic matrix 
sheets spaced by approximately 0.4 microns are created. Crystalline aragonite tablets 
then nucleate on the first layer at uncorrelated random locations and grow while confined 
by the organic matrix sheets. The crystals in this layer continue to grow and fill out 
the space in the layer [1,24,25]. Each nacre tablet in a given layer nucleates and then 
grows until it reaches confluence with a neighboring tablets, so that the resulting tablet 
pattern resembles a Voronoi construction [2,5], with tablet in-plane width of order 5 
microns. The crystal orientation of each tablet is highly probable to be the same as that 
of the tablet directly below its nucleation site [10,11,26], which reflects the presence 
of pores in the organic matrix (typically with diameter ~ 5 — 50 nm) through which 
aragonite crystals can grow. As discussed in [2,8], there is one nucleation site per tablet, 
with an identifiable structure in the organic matrix that is large enough to have one or 
more pores going through it. 

The model that we examine here assumes that growth in a given layer is completed 
before tablets in the succeeding layer are nucleated, and that the positions of nucleation 
sites in all successive layers are nucleated at locations that are uncorrelated with those in 
the preceding layers. The growth rates in the first layer are chosen uniformly at random 
in the interval [1 — 5/2, 1 + 5/2], and the tablets in each layer grow to confluence. With 
probability 1 — e a tablet has the same growth rate as the tablet below its nucleation 
site, while with probability e the tablet is assigned a growth rate chosen uniformly at 
random from the range [1 — 5/2, 1 + 5/2]. Appendix A presents the arguments that 
these choices of the probability distributions are natural, given the geometry of the 
experimental system. The version of the model used here and in Ref. [15] is simpler 
than that examined in Ref. [8], which explicitly models columnar nacre by assuming 
that growth in a given layer is not completed before tablets in the succeeding layer are 
nucleated, so that nucleation sites in successive layers must be within a certain in-plane 
distance and are thus correlated; it is theoretically more attractive because it requires 
the specification of one fewer parameter. Ref. [15] presents numerical simulations of 
the model and experimental data on the pattern of tablet orientations in red abalone 
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Figure 1. Schematic of growth of sheet nacre, which exhibits layer- by-layer growth, 
(a) Scaffolding of organic matrix sheets is created, (b) Crystalline aragonite tablets are 
nucleated at random positions in the first layer of pre-existing scaffolding of organic 
matrix sheets, (c-e) Aragonite tablets are confined by organic matrix sheet in the 
next layer, but grow within the layer until they reach confluence, (f) Aragonite 
tablets are nucleated at random positions in the next layer. With high probability, 
the nucleated tablet has the same crystal orientation as that of the tablet directly 
below the nucleation site [10,11,26]. 

nacre, and demonstrates that numerical simulations of this model using parameter values 
e = 0.015 and 5 = 0.25 yields good agreement with the experimental measurements on 
nacre from red abalone, Haliotis rufescens. 

The next section presents an analytic treatment of the behavior of the probability 
distribution function governing the density of tablets of different orientations as a 
function of the distance from the prismatic boundary. 

2.1. Growth model used for simulations 

This paper reports simulation of a model in which nucleation sites are placed randomly 
on a square calculational domain with open boundary conditions. In each layer, the rate 
of growth of each tablet is chosen at random from a uniform distribution in the range 
[1 — 5/2, 1 + 5/2}. The tablets nucleate simultaneously and grow to confluence before 
successive layers are nucleated. The resulting pattern of tablet boundaries in each layer 
is a multiplicatively weighted Voronoi construction [27-29]. Variations in the in-plane 
extent of each tablet within each individual layer are ignored, so that the volume of a 
given tablet is the product of its area and the layer separation. Because faster-growing 
tablets tend to take up a greater fraction of the area, and because a tablet will nucleate 
with the same velocity as the tablet below it with probability 1 — e, when e is small, 
there is a tendency for faster-growing tablets to take up more and more of the total 
area (and hence volume) as the growth proceeds. The ordering is not perfect, however, 
because some misoriented tablets are nucleated. 
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3. Effective theory for the evolution of probability distribution of tablet 
orientations 

The theory in this section is formulated in terms of crystal orientations, as opposed 
to the formulation in terms of growth rates in the previous section. It is useful to 
consider this formulation because the X-PEEM experiments measure crystal orientations 
as opposed to growth rates, so comparison of theoretical predictions to experimental 
results is facilitated. Motivated by the ordering of tablet c-axis orientations observed 
experimentally, our model posits that variability in the orientations of the c-axes of 
the tablets gives rise to variability in the tablet growth rates, with tablets with c-axes 
aligned perpendicular to the layer plane growing the fastest. 

This section presents and analyzes dynamical equations that govern the evolution 
of the probability distribution describing the number of tablets with different crystal 
orientations. These equations can be solved analytically in the limit that the fraction 
of misoriented tablets is very small, which is the parameter regime relevant to the 
experiments of Ref. [8]. We find that the geometry of the dynamical model leads to 
the result that the distribution of tablet orientations (or, equivalently, growth rates) 
has an extremely sharp peak whose width is exponentially small in e, the fraction of 
misoriented tablets. 

The analytic theory presented here for the evolution of tablet orientations in nacre 
is closely related to mutation-selection models studied in the context of population 
biology [17-22]. The faster growth of tablets of a particular orientation is analogous to 
the process of selection in population biology, where species with higher fitness reproduce 
faster than species with lower fitness. In the growth model examined in this paper, 
after all the tablets in a given layer grow to confluence, the next layer nucleates. When 
a tablet in the next layer nucleates, one of two things can happen. The first, more 
likely, possibility is that the nucleated tablet has the same orientation as the tablet 
directly below its nucleation site (analogous to inheritance, where the descendant has 
the same fitness as the parent), while the second is that the tablet has a randomly chosen 
orientation (analogous to mutation, where the fitness of the descendant differs from that 
of the ancestor by a random amount). Selection in the growth model is reflected in the 
tendency for a larger fraction of the area to be filled with tablets with higher growth 
rates, which increases the probability that a given randomly located nucleation point is 
located over a tablet with a higher growth rate. 

The analytic theory presented here is similar to Ginzburg-Landau theories [30] often 
used in statistical physics [31] to describe phase transitions in condensed matter systems, 
in that the functional forms follow from symmetry and dimensionality considerations 
and may involve unknown coefficients. Even if the numerical values of these coefficients 
are not known, the analytic theory is very useful for obtaining insight into the interplay 
between the various parameters in the problem. However, additional insight is obtained 
if the Ginzburg-Landau parameters can be related to the parameters of the growth 
model described above in Sec. I2.1[ which is done in Sec. H] using a simple mean-field 
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theory and also by comparing the predictions of the analytical model to the results of 
numerical simulations. 

We define ^(7)0^7 to be the fraction of tablets in layer I that are misoriented from 
the layer normal by angles that are between 7 and 7 + d'y. We assume that tablets 
with 7 = have the largest rate of growth, so their share of the area in layer t + 1 will 
tend be greater than in layer £, thus leading to a larger fraction of tablets with crystal 
axes oriented parallel to the layer normal; the function ^(7) governs this tendency. We 
define Xi{l)dl to be the fraction of the area in layer I after its growth is completed that 
is oriented in the range of angles between 7 and 7 + d'y from the layer normal, and that 

Xt(l) = tt^(7)<M7) , (1) 
l\e 

where, for each £, the normalization factor N# is determined by the normalization of the 
probability: 



dl Xi(l) (2) 
N e = I d 1 w( 7 )<M7) • (3) 



The maximum growth rate is at 7 = 0, and one expects w{j) to have a quadratic 
maximum at 7 = 0. We will scale ^(7) so that w(0) = 1. We assume that ^(7) is a 
function only of 7 (in other words, that it depends only on the degree of misorientation of 
the c-axis and not on the orientation of the a and b axes). One expects the dependence 
of w (7) on 7 near its maximum at 7 = to be quadratic, which, if the distribution is 
not too broad, can be approximated Gaussian, 

w(l) oc (l - ^pj » e^ 2 / 2 , (4) 

where a is a numerical coefficient. The tendency for fast-growing tablets to take up 
an increasing fraction of the total tablet area is analogous to the effects of selection in 
population biology, where organisms with higher fitness tend to comprise an increasing 
fraction of the population. 

We then assume that most of the tablets that nucleate in layers above the first layer 
have the same crystal orientations as those of the tablets just below their nucleation 
sites, but that there is a small probability e that a tablet nucleates with a value of 7 that 
is chosen at random from a normalized probability distribution f{y). The nucleation of 
misoriented tablets in the nacre growth model is analogous to the effects of mutation in 
a population genetics model. 

The combination of these growth and nucleation terms leads to a dynamical 
equation governing the behavior of ^(7): 

Ml) = e/(7) + (1 - z)xt(i) 

= e/(7) + ^(l-£)M7W7) • (5) 

To complete the definition of the model, one must specify an appropriate form for 
7(7), the function describing the distribution of angles of misoriented tablet orientations. 
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We will show below that the behavior depends only on the properties of fi^y) as 7 — > 0. 
In Appendix A it is found that the generic behavior for 7(7) for small 7 is for 7(7) to be 
proportional to 7 as 7 — > 0. This result is intuitively reasonable for this system geometry 
because the angular area between 7 and 7 + dj is proportional to 7 as 7 — > 0. At some 
points in the analytic treatment below, we will choose the specific, mathematically 
convenient form 7(7) = fi^e'^ I 2 . This choice does not affect any of the results, 
because the behavior of 7 (7) as 7 — > determines the asymptotic behavior. 

We characterize the behavior of this model using the methods of Refs. [16-18]. 
First we note that when e = 0, so that the co-orientation of tablets in successive layers 
is perfect, this model is easily solved for any initial distribution, ^1(7). Since increasing 
i by one multiplies the un- normalized ^(7) by w(j), it follows immediately that 



(Note that it is sufficient to compute the un-normalized distribution since the 
normalization factor for any given I can always be obtained via Eq. (j2J).) The long- 
time behavior for any w^j) with a single quadratic maximum at 7 = depends only 
on the curvature in ^(7) near 7 = 0. For the specific choice ^(7) = exp(— a7 2 /2), one 
finds <pi(j) oc 01(7) exp(— a(£ — l)7 2 /2). The width of the distribution decreases as the 
square root of the number of layers and becomes arbitrarily narrow as the number of 
layers tends to infinity. 

Now consider the effects of nonzero but small e, so that a nonzero fraction of tablets 
nucleate that are misoriented. The intuitive picture of the process in this regime is that 
the distribution of angles gets narrower unless nucleation of misoriented tablets occurs. 
The misorientations prevent the peak from narrowing indefinitely, so after many layers 
4>e{l) approaches a stationary distribution that does not change as £ increases further. 
Since a fraction e of tablets is misoriented at each layer, one might expect that the peak 
in the distribution narrows for ~ (1/e) layers, so that this naive argument leads to the 
expectation that the eventual width of the distribution should be proportional to e 1//2 . 
However, it is shown below that the subtle interplay between the effects of mutation 
and selection results in a peak width that can be exponentially small in e as e — > 0. The 
high degree of tablet orientation obtained using this mechanism could be a significant 
advantage in the biological context. 

3.1. Steady state behavior of the model 

First we find the fixed point behavior for this model, in which the probability distribution 
function ^(7) approaches a limit (f)*(j) that is independent of £. We expect that this 
fixed point distribution is reached in the limit I — > 00. 

The equation determining the fixed point probability distribution <f>*{pf) is 



Ml) K 4>x (7) (w(t)) 



(6) 




(7) 



with 




(8) 
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The solution to this equation is 



1 " 1^(7) 



with 



.V : / r/- ( /K),r(~)^ 
'0 m=0 



1 - e 

N* 



w(i) 



We now define 



1 -e 

N* 



d 7 (w;(7)) m+ 7(7) , 



and rewrite Eq. ([TO]) as 
1-e 



m=0 



v m I 



(9) 
(10) 

(11) 
(12) 

(13) 



A useful explicit form for the area distributions can be obtained in the parameter 
regime in which e, the fraction of misaligned tablets, is small. When e <C 1, the right 
hand side of Eq. ffl3|) can be of order unity only if v m I m decays slowly for large arguments. 
Because of its definition, v > 0, and it will be seen below that normalization of the 
probability implies that v < 1. Therefore, when m ^> 1, because ^(7) has a single 
maximum at 7 = 0, the integrand in I m from Eq. (TT2T) is very sharply peaked near 
7 = 0. Therefore, the integration interval can be extended to [0, 00], one can assume 
that 10(7) is a Gaussian, 1^(7) = exp(— a7 2 /2), and only the behavior of /( 7 ) for small 7 
is relevant. If the orientations of the misoriented tablets are chosen uniformly at random 
in three-dimensional space, then, as discussed in Appendix A[ ^(7) is proportional to 7 
as 7 — > 0, and for small argument f(j) = fij (with (3 a constant of order unity), so that 
Im — P Jo°° dj 7 exp(— a(m + l)7 2 /2) = (5/ {aim + 1)). Therefore, when e <C 1, 

1 

= f 

a 



v 



It 

n m + 1 

m=0 



(14) 



(Eq. [141 makes it particularly clear that a solution is possible only if v < 1.) Using the 
identity [32] 



E 

k=l 



one finds 



so that 



x 
~~k 

e a 







ln(l — x) 



-ln(l - v) 



N* 





a(l-e)" 


('+■ 


[ fa \ 



(15) 



(16) 



(17) 
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Therefore, 0*(7) can be written 

0*(7) = — 



e/(7) 



1 — ( 1 — exp 



<»(i-«0 

/3e 



exp[— «7 2 /2] 



The width of the probability distribution <p*('j) in 7, estimated by finding the value of 
7 at which 0*( 7 ) = 0*(O)/2, is 



71/2 



2 e _(l_e) Q /(2/3e) 



(19) 



which is extremely small when e is small. 



3.2. Dynamics of the analytic model 



The dynamics of the model defined in Eq. (jSj) can of course be obtained numerically. 
Analytic insight can also be obtained, following Ref. [18], by writing Eq. ((Tj) for the 
fixed point function as 



<t>*{l) 



00 (\ — 

m=0 V 



■w(i) 



(20) 



and comparing this expression to the solution to Eq. 

-1 m— 1 

Mm) = e f ( ci) 



which can be written 



w(7) 0i(7) 



(21) 



En (^»<7) yui 1 - 

m=l n=l x 7 n=l 

If one assumes that the term proportional to ^1(7) is negligible and that the 
normalizations N n can be approximated as being the same, N n = Af, independent 
of n, then one obtains an expression for ^(7): 



<t>t{l) 



m=l ^ 



(22) 



This expression is the sum of contributions that can be interpreted as describing the 
contribution of the population that has undergone a given number of selection events 
since the last mutation, and Ref. [18] shows that it agrees well with numerical solutions 
of Eq. (E]). 



4. Relating the parameters in the effective theory of Sec. [3] to the 
parameters in the growth model of Subsec. 12.11 

Sec. |3] presents an analysis of a theory in which one writes an equation for ^(7), the 
fraction of the area in layer i in which the tablet orientation angle is 7. The model 
analyzed there contains three parameters. The first, a, defined in Eq. HI governs the 
degree to which tablets with small values of 7 grow faster than tablets with larger values 
of 7. The second, e, is the fraction of nucleation sites that have tablets with a random 
orientation instead of the same orientation as the tablet below. The third parameter, 
(3, specifies the width of the distribution governing the distribution of angles of the 
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misoriented tablets. In contrast, the original growth model described in Subsec. 12.11 
has two parameters, e, the fraction of misoriented tablets, and 5, which governs the 
range of in-plane growth speeds of the misoriented tablets. This section describes the 
relationships between the two sets of parameters. 

The value of e is the same in the analytic theory for orientations as in the growth 
model (hence the use of the same symbol). The parameters (3 and 5 are related in a 
straightforward fashion, as discussed below in subsection 14.21 Most of this section will 
focus on the third parameter in the effective theory, a, which we show depends on the 
geometry of the nucleation sites in the growth model in subsection 14.31 We will present 
a simple mean-field theory for estimating the value of a for uncorrelated and random 
nucleation site locations which yields the correct order of magnitude for the value in 
subsection 14.3.11 We compare the mean field predictions with the results of simulations 
of the growth model in subsection 14.3.21 



4-1. Analytic theory for model formulated in terms of tablet growth rates 

In the growth model that is simulated numerically, the tablet orientation angles are not 
considered explicitly because the calculation is formulated using tablet growth rates, 
and the simulation is performed by choosing an initial configuration of nucleation sites 
with a distribution of growth rates, and then allowing the tablets to grow from these 
nucleation sites until they reach confluence. We choose to define a new variable x that 
ranges between and 1, and define Ve(x), the probability distribution describing the 
relative frequencies of tablets in layer i with different values of x [16]. It is natural to 
interpret x as v/v max , where v is the in-plane tablet growth speed and v max is the tablet 
growth velocity at the orientation where this velocity is maximum. 



As discussed above and in Appendix A when the model is formulated in terms of 



orientation angles, the "fitness" function specifying the changes in the fractions of the 
area covered by tablets with different angles of misorientation between successive layers 
is expected to have a quadratic maximum at 7 = 0. Because the tablet growth velocity 
itself depends quadratically on 7 near 7 = (again, because misorientations by 7 and 
—7 are equivalent and so yield the same growth velocity), the "fitness" function that 
specifies the change in relative area of the different values of x, w(x), depends linearly on 
x near x — 1. Recalling that the probability distributions in the model are normalized, 
so that the overall scale of w(x) is arbitrary, and that the behavior is dominated by 
the behavior near x = 1, two ways of parameterizing this dependence are to (1) fix the 
value of w(l) = 1 and specify the slope w', or (2) to write w(x) = x^ and specify £. 
These two forms are equivalent near x — 1, with w'(x) = £; we will choose to use the 
power-law form w(x) = x^. 

To obtain additional insight into the relationship between the two formulations of 
the model, we reformulate the analytic theory in terms of the variable x (this formulation 
is very similar to that in Ref. [16]). Assuming that tablets with normalized growth rate 
x grow to have an area that is proportional to x^ for some £, and assuming that there 
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the probability of nucleation of a misoriented tablet is e, the probability distribution for 
the tablet growth velocities in layer £, Ve, obeys 

P e+ i(x) = eg{x) + (1 - e){l/N t )a*V t {x) , (23) 

where g(x) is the probability distribution for the misoriented tablets and J\fi = 
J* dxx^Ve(x). Using methods precisely analogous to those in the previous section, we 
write the equation for fixed point reached at large £, V*(x): 

eg(x) 



V*(x) 



1 _ (lz£l^ 
1 Af* 



v / '\ el x "' 



m=0 



where W* obeys 



»i 

M*= I dxV*(x) 





dxeg(x)J2[^J F 1 ^) • (25) 

When e is small, then the integrals are dominated by the region where x is close to 
unity, and the x-dependence of g(x) can be neglected, yielding 

m=0 x ' 

When e is small, so that v is close to 1, this sum is dominated by the region in which 
£m ^> 1, and, using Y^k=i xk /k = ~ l n (l — x )-> one obtains 

M*= -^Infl- 1 —^] . (27) 



Recalling that Af* is extremely close to 1 — e, one obtains 

l-^-^taCl-i^, (28) 

so that 



£ V AT* 



A/~* = (1-6) ( 1 + exp- ( vl ^ ) ) (29) 



and 



V*(x) 



eg(x) 



^1 — exp 



(i-«Qg 

69(1) 



(30) 



This form is consistent with Eq. ffl8l) . Comparing the forms, one can see that the 
parameter a in the effective model enters into Eq. ffTSl) in a way similar to that in which 
£ enters into Eq. (130]) . and that (3 enters into Eq. (fl8|) in a way similar to that in which 
g{l) enters into Eq. (|30|) . 



Model of crystal orientation in nacre 



12 



4-2. Relationship between parameters (3 and 8. 

The parameter (3 in the analytic model and the parameter 8 in the growth model 
each describe the probability distribution governing the misoriented tablets. As seen 
above, because of the action of the "selection" in the model, the key property of these 
distributions is their behavior near 7 = and x — 1. Therefore, we work out the 
relationship between the two formulations of the model in this region. 

For the growth model, the distribution of the in-plane growth speed v is chosen to be 
uniform in the range [1 — 5/2, 1+5/2], so the probability density in terms of the variable v 
is 1/8. To convert to x, the normalized growth speed, we note that the maximum growth 
velocity is 1+5/2, so that the probability density for x is (l/8)(dv/dx) = [1+5/2)/ 8. For 
the analytic formulation of the theory, choosing this distribution to be f(j) = f3^e~ 131 ' 1 ^' 2 
(the coefficient is determined by normalizing /(t)^7 — 1, and assuming that (3 is 
large enough that setting the upper limit of the integration at 00 introduces negligible 
error), the probability density near 7 = is f(j) — > ^7. Near 7 = 0, the variables 7 
and x are related by x = 1 — Aj 2 , for some constant A, so near 7 = we can write 

^7/(7) ~ dx(l + 8/2)8 

dx(dj/dx)(@j) « dx(l + 8/2)/ 8 

((3 1 )/(2A 1 ) * {1 + 8/2)/ 8, (31) 

leading to the identification (3/ [2 A) = (1 + 8/2)/ 8. This relationship between the 
probability distributions of the orientation angles and the growth speeds involves not 
just (3 and 8 but also a new constant A, that specifies the relationship between the 
change in orientation angle and the change in growth speed. However, this constant A 
also enters into the "fitness" or "selection" term that is discussed in the next subsection 
in precisely analogous fashion, so that the relationship between ratios a/f3 and £/g(l) 
in Eqs. (TIB]) and Eq. ( 1301) does not depend on A. 

4-3. Dependence of a on geometry 

This subsection discusses how the parameter a in Eq. H] is related to the growth model 
as described in Subsec 12. 1[ It can be seen that a does not depend on e and 8 by noting 
that in the limit of e — > 0, when no misoriented tablets are nucleated, both e and 8 drop 
out of the model altogether, while the selection term in which a appears still operates. 
We will see that the value of a is essentially determined by geometrical considerations 
only. 

The value of a in the analytic model depends on the geometric arrangement of 
nucleation sites in the growth model. This is because the areas of the tablets after 
growth is completed depend on the differences between growth rates of neighboring 
tablets as opposed to depending on individual growth rates themselves; two neighboring 
tablets growing at the same speed will not lose area to each other, no matter what that 
rate happens to be. Changes in the fractions of area covered by tablets with different 
growth rates arise only when fast-growing tablets are next to slow-growing tablets. In 
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fact, if one considers a one-dimensional case of iV tablets in which the growth rate of 
tablet j is g + cj, where g + and c are constants, it is easy to see that there is a set of 
locations for the nucleation sites such that the size of every tablet stays the same except 
for the tablets at j = 1 and j = N. 

Here we consider the specific case of a collection of uncorrelated, randomly chosen 
nucleation sites. For this situation, neighboring tablets are likely to have different growth 
speeds, so that the differential growth rates will lead to fast-growing tablets increasing 
their share of the area. We first present a simple mean field estimate for the resulting 
value of a, and then present numerical computations that can be used to extract its 
value. 

4-3.1. Mean field theory for a for random nucleation sites In this subsection we 
construct a simple mean-field theory for estimating a. In the next subsection we 
will compare the mean-field theory results with the results obtained using numerical 
simulations. 

In subsection 14.11 we showed that without loss of generality, one can choose the 
form for the selection term in the growth model for the relative area covered by tablets 
with normalized growth speed x to be x^. Here we present a simple mean-field theory 
for calculating £ using arguments very similar to those typically used in statistical 
physics [33]: We assume that when a tablet grows, it does so in the presence of neighbors 
with mean properties, and take account explicitly of the fact that faster-growing tablets 
will tend to reach confluence with their neighbors sooner than slower-growing ones. One 
assumes that when a given tablet grows, it does so in the presence of neighbors that 
have average properties. 

Given that a tablet has normalized growth speed x = v/v max and that the mean 
normalized growth speed is x = v/v max , then we wish to calculate the area covered by 
the tablet after it has grown to confluence with its neighbors. For simplicity, we will 
consider nucleation sites that are placed on a regular lattice (or alternatively, consider 
the scaling of the behavior without worrying about numerical coefficients). For two 
nucleation sites separated by d with normalized growth rates X\ and x, the time for the 
tablets to grow to confluence is proportional to d/(x\ +x), so the area of the tablet in 
question should scale as X\ 2 J(x\ + x) 2 . Therefore, the ratio of the area taken up by the 
tablet we are considering to that of the "mean-field" tablet (which we assume has area 
that scales as x 2 / {x\ +x) 2 ) is X\ 2 jx 2 . Approximating x as constant (which is reasonable 
because the probability distribution for x is asymmetric, monotonically decreasing from 
x — 1, with a mean very near 1) yields the result that the area growth for a normalized 
velocity X\ is proportional to xf, so that £ = 2. 

A modification of this mean field theory, that yields a value for £ that agrees better 
with the results of the numerical simulations below, is obtained by noting that while the 
area of the tablet under consideration is proportional to xi 2 /(xi + x) 2 , the area of the 
mean-field neighbors varies little when x\ changes because their interactions with their 
other neighbors are not affected. Assuming that the change in areas of the "mean-field" 
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tablets can be neglected altogether yields the value £ = 1. This result can be seen in 
several ways: first, the mean- field argument given just above assumes that increasing 
x\ leads to a decrease in the "mean-field" tablet area that scales the same way as the 
increase in the area of the "non- mean- field" tablet, so ignoring the former contribution 
leads to a result that is the square root of the one obtained previously. Alternatively, 
one can argue as follows: since the area growth of the "non-mean-field" tablet growing 
with normalized velocity X\ is proportional to x\/(x\ + x) 2 , we can write x± — 1 — p±, 
x — 1 — p, assume that p\ and p are each <C 1, and fix the proportionality constant 
by requiring the relative area to be unity when x± — 1. When p is close to unity, this 
procedure also yields the result £ = 1. 



4-3.2. Numerical evaluation of parameters in effective theory for synchronous nucleation 
and random nucleation geometry For the simulations presented in this section, the 
initial distribution of tablet growth velocities tablets are chosen to be uniform in the 
range [1 — 5/2, 1 + 5/2], or, equivalently, g(x) = 1/5 for x in the range [(1 — (5/2) / (1 + 
5/2), 1]. Because the focus is on determining the value of a in the selection term of the 
effective model, the parameter e is set to zero in the simulations in this section. 

The comparison between mean field theory and simulation is done by comparing 
x i/2 {(■), the median growth speed found in the simulation at layer £, to the median of 
the distribution found using the analytic mean-field theory, which satisfies 

V e (x) = 1/2 . (32) 

'o 

When e = 0, the dynamical equation for Ve(x) is 

V e+1 (x) (xx^V e (x) , (33) 

which has the solution 

V t (x) = x^V (x) . (34) 

When £ is large, the x-dependence of Vo(x) can be neglected, and one obtains Vifx) oc 
a£ , yielding a normalized distribution Ve{x) = + l)x^ . The median value of this 
distribution, Xi/ 2 (£), is 



f 

Jo 



x 1/2 (e) 



1 \ e<+i 
2 



ln2 . . 

(35) 



^£+1 

Therefore, when £ is large, the slope of the plot of the median value X\/2{£) versus 1/ £ 
approaches ln(2)/£, providing a method for extracting the value of £ from numerical 
simulations. 

The procedure used for the numerical simulations is very similar to that used 
for those presented in Ref. [8], but specialized to the limit in which successive layers 
nucleate much more slowly than individual layers grow to completion. Each layer in 
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the simulation is a square of dimension 400 x 400 (the length unit is arbitrary, and 
is adjusted to agree with measured tablet widths when numerical results are compared 
with experimental data in Ref. [8,15]) with nucleation sites whose sites are independently 
and randomly chosen from a probability distribution with uniform spatial density. Open 
boundary conditions were employed. The growth rates in the initial layer are chosen 
uniformly at random in the interval [1 — 5/2, 1 + 5/2] (resulting in normalized growth 
rates in the interval [(1 — 5/2) /(l + 5/2), 1]). In these simulations, the parameter e is 
set to zero (no misoriented tablets are nucleated). 

4-3.3. Numerical results As discussed in Sec. I4.3.2[ the value of £ that determines 
the strength of the "selection" can be extracted by computing the median value of the 
in-plane tablet growth velocity, x me d(£) as a function of layer number £, when e = 0. 
Eq. (|35|) predicts that the quantity 1/(1 — x me d(£)) depends linearly on £, and the slope 
of this dependence is £/ln2. For the uniform distribution of growth rates used in the 
numerical simulations, the initial value of this quantity, 1/(1 — x me d(l)), is 1 + 2/5. 

Figure 2 is a plot of the quantity 1/(1 — x me d(£)) as a function of layer number £ 
for different values of 5. The graph also shows the dependence expected for 5 = 0.3 
using the two different mean-field results obtained above, £ = 2 and £ = 1. There is 
substantial scatter in the numerical data, but it appears that the behavior is consistent 
with the value £ = 1 but not with the value £ = 2. 
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Figure 2. Plots of 1/(1 — x me d(t)) versus layer number £, obtained by numerical 
simulation of the model, where x me d{£) is the median value of normalized growth 
speed x in layer £. The number of nucleation points per layer and the value of S 
for each run are shown in the legend. For the uniform distribution of growth rates 
used in the numerical simulations, the initial value of the ordinate, 1/(1 — x me d(l)), is 
1 + 2/(5. The mean field theories described in the text yield a linear dependence for this 
quantity, one with slope 2/ In 2 « 2.88, corresponding to the value £ = 2 (a dotted line 
with this slope is shown on the plot), and the other with slope 1/ In 2 w 1.44 (dashed 
line on the plot). The numerical results exhibit considerable scatter, but appear to be 
consistent with the value of £ = 1 and not the value of £ = 2. 
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5. Discussion 

In this paper we have detailed investigations of a model for the development of 
orientational order of the aragonite tablets in nacre. The model is investigated both 
by direct numerical simulation and by analysis of the evolution of the probability 
distribution governing tablet orientations that is analogous to mutation-selection 
theories studied in population biology. The relationships between the parameters in 
the version of the model used for the numerical simulations and those in the version 
used for the analytic work are investigated. 

The model simulated in this paper can be extended in many ways to model nacre 
growth more realistically. In Ref. [8] the model is modified to take the columnar growth 
of abalone nacre into account, by restricting each nucleation site to be within a given 
distance from a nucleation site in the layer below, as opposed to the model used in this 
paper, in which the nucleation sites in each layer are chosen independently and randomly 
in the simulational domain. Compared to the version in Ref. [8], the model examined 
here has the advantage that it has one less parameter. Because the analytic theory is 
constructed using very general symmetry considerations, it is clear that the qualitative 
behavior of a model implementing columnar growth is the same as for the model studied 
here. Similarly, the growth model can be generalized in other ways (for instance, by 
implementing more realistic microscopic growth rates with explicit modeling of three- 
dimensional growth or by considering asynchronous nucleation) will still yield the same 
analytic model, though with possibly different values for the model parameters. 

The relationships between the parameters in the analytic theory and in the growth 
model presented in Subsec. 12.11 were investigated using a mean field theory and by 
numerical simulation. Two mean field theories were constructed to estimate the 
parameter describing the strength of the tendency for faster-growing tablets to take 
up more area than slower-growing ones. Numerical simulations of the growth model 
yield a result for this parameter that agree reasonably well with one of the mean field 
theories. Development of a more sophisticated analytic theory for calculating the value 
for this parameter is an interesting open problem. 

6. Conclusions 

This paper investigates theoretically a model for the development of orientational 
order of the aragonite tablets in nacre, or mother-of-pearl. Motivated by experiments 
demonstrating the this ordering develops over many tens of layers, the model assumes 
that the tablet growth rates depend on the orientation of their c-axes, with tablets with 
c-axis orientation normal to the layer plane growing the fastest. This model is closely 
related to mutation- select ion models used in population biology. 

A combination of analytic and numerical results were applied to analyze the model. 
It is shown that the model yields a degree of ordering that is extremely good when the 
parameter e, that governs the fraction of tablets with nucleate with misoriented c-axes, 
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is small. 
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Appendix A. Form of functions used in analytic model 

This appendix discusses the forms of the functions w(j) and f(j) that are used in 
Sec. [3] to specify the model for the evolution of the probability distribution of tablet 
orientations. The function w( , y) specifies the tendency of the area covered by faster- 
growing tablets to increase, while 7(7) specifies the probability distribution function for 
the angles of the misoriented tablets. We assume that the growth rate is isotropic in 
the a-b plane, so that the tablet growth velocity depends only on 7, the angle between 
the c-axis orientation of the tablet and the layer normal. 

A 




Figure Al. Schematic of geometry, illustrating that assuming that the c-axes of 
misoriented tablets are equally likely to have any orientation, then the function / (7) 
introduced in Sec. [3] is proportional to 7 as 7 — > 0. 

For the function ^(7), our assumption that the growth speed has no azimuthal 
dependence implies that the growth speed of a tablet with crystal axis misoriented from 
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the layer normal by an angle 7 is a function of cos(7) and has a quadratic maximum 
at 7 = 0. When the parameter e governing the fraction of misoriented tablets is small, 
the dynamics is controlled by the behavior for small 7. As discussed in the main text, 
it is often convenient to choose 117(7) to be a Gaussian, which makes calculation of its 
dynamical effects particularly simple. 

Now we discuss the function 7(7) that describes the distribution of misoriented 
tablets. If one makes the simplest assumption that the orientations of these misoriented 
tablets take on each angle with equal probability, then the probability distribution for 
7, the angle between the tablet c-axis orientations and a fixed external axis (here, the 
direction of the layer normal) is proportional to sin(7), since the area of a spherical 
surface between angles 7 and 7 + dj is proportional to sin(7)rf7, or dcos(7). The 
relevant geometry is shown in Fig. 3. 

As discussed in the main text, the model can be formulated entirely in terms of 
tablet growth rates, and we use such a form in our numerical simulations. We define 
the variable x to be the normalized tablet growth rate, with x — 1 the maximum value. 
Because (1) we assume that the tablet growth rate exhibits a quadratic maximum at 
7 = 0, and (2) the area between the angles 7 and 7 + dj is proportional to 7 as 7 — > 0, 
in terms of x, the growth rate is maximum at x — 1 and depends linearly on x near 
x — 1. Similarly, the function describing the relative change in area taken up by tablets 
with c-axis orientation 7, 117(7), has a quadratic maximum as a function of 7, so the 
analogous function expressed in terms of the growth rate will depend linearly on x near 
x = l. 
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